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Abstract. We study spin 3/2 fermionic cold atoms with attractive interactions confined in a one-dimensional optical 
lattice. Using numerical techniques, we determine the phase diagram for a generic density. For the chosen parameters, 
one-particle excitations are gapped and the phase diagram is separated into two regions: one where the two-particle 
excitation gap is zero, and one where it is finite. In the first region, the two-body pairing fluctuations (BCS) compete 
with the density ones. In the other one, a molecular superfluid (MS) phase, in which bound-states of four particles form, 
competes with the density fluctuations. The properties of the transition line between these two regions is studied through 
the behavior of the entanglement entropy. The physical features of the various phases, comprising leading correlations, 
Friedel oscillations, and excitation spectra, are presented. To make the connection with experiments, the effect of a 
harmonic trap is taken into account. In particular, we emphasize the conditions under which the appealing MS phase 
can be realized, and how the phases could be probed by using the density profiles and the associated structure factor. 
Lastly, the consequences on the flux quantization of the different nature of the pairing in the BCS and MS phases are 
studied in a situation where the condensate is in a ring geometry. 

PACS. 03.75.Mn Multicomponent condensates; spinor condensates - 71.10.Pm Fermions in reduced dimensions 
(anyons, composite fermions, Luttinger liquid, etc.) - 71.10.Fd Lattice fermion models (Hubbard model, etc.) 

1 Introduction and magnetism, which stems from the presence of the different 

internal states, can be investigated. Experimentally, three 
Recent experimental progress achieved in trapped ultracold component Fermi gases can be created by trapping the three 
atomic gases provides a great opportunity for exploring the lowest hyperfine states of ^Li atoms in a magnetic field, or 
physics of strong correlations in clean systems, thanks to the ^Y considering "K atoms. In addition, the magnetic field 
tunability of interactions using optical lattices and Feshbach dependence of the three scattering lengths of ^Li is known 
'<^ ■ resonance. A large number of interesting phenomena of con- experimentally and can be tuned via Feshbach resonance l32l 
densed matter physics and nuclear physics is then expected to which opens for the experimental realization of a three- 
be accessible in the context of ultracold atomic gases HI. A component fermionic lattice model. In fact, such a degenerate 
prominent example is the observation of the Mott insulator- Fermi gas has been reaUzed experimentally very recently [33J, 
superfluid quantum phase transition with cold bosonic atoms in and a four-component Fermi gas could also be achieved using 
an optical lattice |2|, and its possible fermionic analogue, the '^ atoms |I34|. 

Mott insulator-metallic phase transition, recently investigated The existence of these internal degrees of freedom is ex- 

in a two-component Fermi gas fSj. A second breakthrough is pected to give rise to some exotic superfluid phases. In this re- 

the trapping of a two-component Fermi gas and the study of the spect, a molecular superfluid (MS) phase might be stabilized 

crossover from fermionic superfluidity of Cooper (BCS) pairs where more than two fermions form a bound state. Such a 

to Bose-Einstein condensation of tightly bounded molecules m state might be relevant to several topics in physics. For in- 

|5]|6l. stance, the quark model of nuclear matter at low density de- 

The superfluid behavior of multicomponent Fermi gases scribes nucleons as three-fermion bound states. Such a trionic 

with more than two hyperfine states might also lead to inter- phase has been found in one-dimensional integrable fermionic 

esting properties that have been explored recently [7.,8,9i.l0, model with three colors |7| and its emergence in the context 

[niiniinninMIMnillSlIllEilll]!^ of three-componem ultracold fermions has been discussed re- 

l29ll30ll3TI ■ In particular, the interplay between superfluidity cently M22il26ll2"7ll29lf30ll . The possibility that superfluidity is 
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sustained by a condensate based on four-fermion bound states 
(quartet) might be also explored in cold atomic physics lITTl 
[14, 15 16,26 31 1. Such a superfluid behavior has already been 
found in very different contexts such as nuclear physics for 
instance, where a four-particle condensate-the a particle-is 
known to be favored over deuteron condensation at low den- 
sities 1 35, 36 1 . Such a quartet condensation can also occur in 
semiconductors with the formation of biexcitons |37|. A quar- 
tetting phase, which stems from the pairing of Cooper pairs, has 
also been found in a model of one-dimensional (ID) Josephson 
junctions |38| and in four-leg Hubbard ladders [39 1. 

In this paper, we will investigate the low-energy properties 
of (hyperfine) spin-3/2 (i.e. four-component) fermionic cold 
atoms confined in a one-dimensional optical lattice in light of 
the possible formation of a quartetting phase. Due to Pauli's 
principle, low-energy s-wave scattering processes of spin 3/2 
fermionic atoms are allowed in the singlet and quintet chan- 
nels, so that the effective Hamiltonian with contact interactions 
reads |[lp,40| 



^=-*E[^a,«^"=^+l+h-C-]-E 



^iUi 



^0 2Z -^OO.i^OO.i + C^2 2^ ^2m,»-f2m,i , 



(1) 



where c*^ ^ is the fermionic creation operator at site i, in one 
of the a — ±1/2, ±3/2 hyperfine states. The on-site density 

operator is denoted by n^ = X)a "^I i*^".*- ^^^ chemical po- 
tential Hi can be uniform (called ^ for grand-canonical Quan- 
tum Monte-Carlo calculations), inhomogeneous in presence of 
the trap, zero for DMRG calculation (canonical ensemble). For 
convenience, the lattice spacing is set to unity. Singlet and quin- 
tet operators in Eq. ([T]i are defined using Clebsch-Gordan coef- 
ficients 

a/3 

For instance, the spin 3/2 on-site singlet operator reads pJq ^ — 

Pi = cI/T c „/„ — cj ,„ .c* -1 ,„ .. A convenient way to rewrite 
the Hamiltonian is to express it in terms of the density and sin- 
glet pairing operators: 



^ = -^E[''L^"''+i+'^-'^-]"E''» 



V 



Y^nUvY^PlP, 



(2) 



with U — 2U2 and V = Uq ~ U2- This model has an exact 
SO(5) symmetry II41II42I . and, for the fine-tuning Uq = U2 (or 
y = 0), a SU(4) symmetry. In the latter case, the Hamilto- 
nian reduces to a Hubbard-like Hamiltonian with only on-site 
density-density interactions. It resembles the usual SU(2) Hub- 
bard model, but with four colors instead of two and we refer to 
it in the following as the SU(4) line. Similarly, we refer to the 
[/ = and V < line as the BCS line since the singlet pairing 
is naturally favored in this regime. The model (|2|i has essen- 
tially three physical parameters: the density of particles n, and 
the two interactions U/t and V/t in units of the hopping t (set 



to one in the following). Experimentally, the interacting param- 
eters can be varied by tuning the scattering lengths (for instance 
to negative values) and the depth of the optical lattice |4|. 

In the homogeneous situation (i.e. in absence of the har- 
monic trap), the phase diagram of model Q at zero tempera- 
ture has been investigated by means of low-energy approaches 
[14,15,311 and numerical calculations |26,43| such as the 
density-matrix renormalization group (DMRG) technique 1441 
I45ll46l and Quantum Monte-Carlo (QMC) simulations II47II48I 
l49l . Away from half-filling, there are two very different spin- 
gapped phases which are separated by an Ising quantum phase 
transition. In the first one, for instance along the SU(4) line 
with U < 0, the BCS singlet-pairing instability is suppressed. 
The leading instability is an atomic-density wave (ADW) with 
wave-vector 2/ci? (kp being the Fermi wave-vector) or a quar- 
tetting one. In particular, at sufficiently low-density, a dominant 
MS instability emerges which marks the onset of the quartet- 
ting phase. In the second spin-gapped phase, basically obtained 
along the BCS line, the 2fci?-ADW instability has now a short- 
range behavior and BCS singlet pairing competes with a molec- 
ular density-wave (MDW) with a Akp wave-vector. 

In this paper, we give more details on the large-scale numer- 
ical calculations which have been used in the short papers 1261 
43] and bring several new results. In this respect, we present the 
phase diagram of model Q in absence of the trap for a generic 
filling which is not one atom per site as in Ref . I*??] . Moreover, 
the Friedel oscillations and excitation spectra are studied, to- 
gether with the quantum phase transition between the two spin- 
gapped phases from the behavior of the entanglement entropy. 
We also introduce a simple observable, the molecules fraction, 
which could be useful for experiments. Then, we investigate the 
inhomogeneous situation and the effet of a harmonic confining 
potential on the quartetting phase in order to make contact with 
future experiments in spinor fermion ultracold gases. Finally, 
the nature of the flux quantization in the BSC and MS phases 
is analyzed in a ring geometry. 

The paper is organized as follows. In Section |2l we recaU 
the main results obtained within the low-energy approach and 
we describe the technical details of the three numerical meth- 
ods used in this work. Section[3]presents our main results con- 
cerning the phase diagram and the physical properties of the 
phases of model (|2|i in the homogeneous situation. The exper- 
imental signatures of the quartetting phase are then discussed 
in Section|4]which includes, in particular, the effect of the trap. 
Finally, our concluding remarks are summarized in Section|5] 



2 Low-energy and numerical approaches 

2.1 Low-energy approach 

In this section, we recall the main results of the low-energy 
approach 1 14, 15,31,501 on the behavior of the different order 
parameters that identify the possible phases of model (|2]). For 
a generic density, the low-energy Hamiltonian separates into 
two commuting pieces: a density and (hyperfine) spin part. 
This result is nothing but the famous "spin-charge" separa- 
tion which is the hallmark of ID incommensurate electronic 
systems M51II52II . The U(l) density fluctuations remain gapless 
while the spin part is fully gapped for parameters with either 
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Phases | 


ADW (K < 2) 
MS {K > 2) 


BCS (K > 1/2) 
MDW (K < 1/2) 


Correlator 


exponent 


wave-vector 


exponent 


wave-vector 


Q{x) 


2/K 





2/K 





P{x) 


exp. 





1/{2K) 





N2kA^) 


K/2 


2kF 


exp. 


2kF 


NikA^) 


2K 


Akp 


2K 


-ikp 



Table 1. The possible phases, obtained by means of the low-energy 
approach, of spin-3/2 cold atoms with attractive interactions; the sym- 
bol exp. denotes a correlation with an exponential decay and the 
other correlations have a power-law behavior; N2kp (x) (respectively 
^^4^;^^ (x)) corresponds to the 2kp (respectively 4kp) part of the den- 
sity correlation function. 



U otV negative. In the terminology of ID electronic systems, 
we have the stabilization of a Luther-Emery liquid phase ifSTl 
|52]| . However, in contrast with the standard F = 1/2 (i.e. two- 
component) fermions, two very different Luther-Emery liquid 
phases emerge here, which are separated by an Ising quantum 
phase transition when U = V < Oin the weak-coupling limit. 
In the first one (dubbed ADW/MS region), a representative 
being the SU(4) line with [/ < 0, the BCS singlet pairing with 
order parameter Pi displays a short-range behavior. The lead- 
ing instability of this region corresponds to the order parame- 
ter which has the slowest (power-law) decaying correlations at 
zero temperature. The natural candidates in this first phase are 
the density correlation A^(x) = (riirii+a;) and the quartet corre- 
lation (3(a;) = {QiQl+J withQi = C3/2,iC_3/2,jCi/2,iC_i/2,,. 

The latter instability is rather natural since when U < and 
V = 0, the density-density term in Eq. (|2|i has a tendency to 
form on-site molecules (the quartets) made of four particles in 
all different hyperfine states. The leading asymptotic behavior 
of these correlations can be computed within the low-energy 
approach and the results are summarized in Table[T] The power- 
law decay of these correlation functions depends only on the 
Luttinger parameter K which is a non-universal function of 
the interaction parameters and the density n. A perturbative es- 
timate, valid in the weak-coupling regime |V|/i, \U\/t ^ 1, 
gives 



K ^ l/y/l + [V + 3U]/7rvF , 



(3) 



with the Fermi velocity vp — 2isin(7r?i/4). Within the same 
approximation, the sound velocity of the gapless density mode 
reads: 

u = VF^l + [V + 3U]/nvF . (4) 

For non-interacting fermions, we have K — 1 and u ~ vp. 
From Eq. (O, we see that attractive interactions increase the 
Luttinger parameter and that an artificial divergence, which oc- 
curs when the denominator vanishes, signals the breakdown of 
the perturbative calculation. From Table [T] we deduce that the 
leading instability is the 2/ci?(= Trn/2) ADW one for K < 2, 
whereas the quartetting MS phase is stabilized when K > 2.ln 
the latter regime, we have an exotic Luther-Emery liquid with 
a confinement of pairs (that would be objects with charge 2e in 
a context of charged particles) and the emergence of quartets 
(similarly, objects with a 4e charge). A related Luther-Emery 



phase has been found in a totally different context correspond- 
ing to the formation of multi-magnon bound-states in the spin- 
1/2 J1-J2 Heisenberg chain under magnetic field |53|. Inside 
the ADW/MS region, there is no sharp quantum phase transi- 
tion and only a smooth crossover. In this respect, it might be 
interesting to observe that there is a continuity between weak 
and strong coupling regimes in this region. Indeed, the higher- 
harmonics in the quartet correlation can be estimated by means 
of the low-energy approach: 

^-{2/K+K/2) 

+ Ccos(4A:pa;)x-2(^+i/^), ^^^ 

A, B, C being non-universal amplitudes. On the other hand, 
along the SU(4) line at small densities and strong attractive 
U, the physics is essentially governed by hard-core bosons 
bi ^ Qi with repulsive interactions. The bosonic correlation 
function of this model is known from the harmonic-fluid ap- 
proach IMlEll : 



Q{x) ^A x-'^l^ + Bcos{2kFx) 






(6) 



where po ~ n/4 is the density of the bosons and Kf, is the 
underlying Luttinger parameter From Eqs. (|5]l and (|6]l, we thus 
observe that there is a continuity between weak and strong cou- 
pling regimes with K — 4Ki,. In particular, we also deduce an 
upper bound for the Luttinger parameter K: -K^max = 4 since 
the value Kb = 1 for non-interacting hard-core bosons 15411551 
is expected in the limit of vanishing densities. 

In the second spin-gapped region (called in the following 
BCS/MDW region), obtained for instance along the BCS line, 
the pairing term in Eq. (|2| stabilizes the order parameter Pi of 
the Cooper pairs. Now, the 2k p ADW instability is a strongly 
fluctuating order since the 2fci? part of the density correlation 
has an exponential decay. As seen in Table[Tl the competing or- 
ders in this phase are the BCS instability with equal-time cor- 
relations P{x) ~ {PiPl^x) and 4kF{= Trn) ADW operator. 
A BCS phase is stabilized for K > 1/2 which is analogue to 
the standard Luther-Emery phase of spin- 1/2 electrons 15111521 . 
For K < 1/2, a MDW phase, which is characterized by a 4kF 
oscillation of the density fluctuations, is predicted to emerge. 
For a generic filling, we expect no quantum phase transition 
between BCS and MDW phases but a smooth crossover For 
the commensurate filling of one atom per site (n = 1), we have 
shown in Ref. |43| that a Mott transition occurs and that the 
MDW phase is replaced by a Mott-insulating phase with bond 
ordering. 

In summary, we observe that the nature of the phases 
found within the low-energy approach are governed by the non- 
universal Luttinger parameter K which is a function of the den- 
sity n and the interactions U/t, V/t. It is thus crucial to have 
a reliable evaluation of K. Since model ^ is not integrable in 
the generic case, numerical calculations of this parameter are 
required. 

2.2 Numerical methods 

We use three different numerical methods to investigate the 
phase diagram of model (|2|i: mainly the density-matrix renor- 
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malization group (DMRG), but also the exact diagonalization 
(ED) and quantum Monte-Carlo (QMC) techniques. 

DMRG calculations were performed at zero temperature 
with open boundary conditions (OBC) using an exact mapping 
of model (|2|i onto a two-leg SU(2) Hubbard ladder model with 
special couplings (here, the spin index can take only two values 

a = ±1/2): 



3 Phase diagram 

This section gathers results on the phase diagram for a generic 
density, i.e. a density for which no commensurability effects 
are expected, and which is sufficiently low to realize the MS 
phase on the SU(4) line ll26l . We choose n = 0.75. We first ex- 
plain how the Luttinger parameter K is computed numerically, 
before giving more details on the physics of each phase. 






(7) 



We use L as the label for the ladder couplings, and _L for cou- 
plings between the two chains and || for couplings along the 
chains. (3 = 1, 2 is the chain index, rii^p^a = c\ ^ ^Ci^p^a, ^ijj 
is the spin operator, and ri; ,g = ^^ ni^p^a the local density. 
For the hoppings, we have A = t and t\^ = 0. For the on-site 

interaction, C/*- — U. For the next-nearest neighbor density- 
density interaction on rungs, V\^ = U+V/2, and also a Heisen- 
berg coupling on the rungs ,/^ — —2V. All other couplings 
are equal to zero. In this mapping, the local number of states 
per site is strongly reduced as it is 2^ = 4 for a SU(2) Hub- 
bard site and 2'* = 16 for a spin-3/2 Hubbard site. Symmetries 
are used to fix the total number of fermions to Nf and the to- 
tal z-component of the spin to zero. We have typically kept 
1000 states of the reduced density matrix, and sometimes up 
to 1400. Convergence depends on the regions of the phase dia- 
gram, and is harder with the trap. On the SU(4) line with large 
\U\, the convergence is very good as the physics is essentially 
the one of hard-core bosons. The discarded weight typically 
ranges from 10~^^-10~^ when both interactions are negative 
(and not to small) to 10^^ if one is positive or small. Moreover, 
the discarded weight decreases with density so that simulations 
become easier and more accurate in this regime. 

On the SU(4) line with total S^- — 0, the numbers of par- 
ticles Ncr per specie are independently conserved. Therefore, 
with an appropriate choice of boundary conditions and N^ (for 
instance periodic boundary conditions and N^ odd), the par- 
ticles do not experience any statistics so that, by means of 
a Jordan- Wigner transformation, the model is strictly equiva- 
lent to a hard-core boson model on a four-leg ladder for which 
chains are only coupled via a density-density interaction term 
V\^ ^ U between all chains. Such a bosonic model has no 
sign problem and can be efficiently simulated by QMC tech- 
niques such as the Stochastic Series Expansion (SSE) algo- 
rithm P47''481. We use the ALPS software implementation of 
the SSE algorithm |56 57|. Note that, contrai'ily to DMRG, the 
algorithm works in the grand-canonical ensemble and at finite 
temperature. Away from the SU(4) line, Fermi statistics can- 
not be avoided and therefore, we have also used a determinan- 
tal QMC algorithm (DQMC), which has no sign problem over 
a relatively wide range of parameters 1.43 J . For this algorithm, 
we have used the projector approach that provides ground-state 
properties with a fixed number of particles B9l . 



3.1 Extracting the Luttinger exponent 

As one can see from Table [T] the Luttinger exponent K can 
be extracted from algebraically decaying correlation functions. 
For instance, the quartet (or MS) correlations Q{x) gives ac- 
cess to 2/K in all regions of the phase diagram. They can be re- 
liably computed with DMRG if the number of state kept is suf- 
ficiently large. As data are computed on finite and open chains, 
there is no translational invariance and all correlators depend 
on both positions of the sites. For instance, Q{x) = {QiQ]^^) 
will actually depend on i. We fix i = m = L/2 to be at the 
middle of the chain and control finite size effects using re- 
sults from conformal theory |55|. By denoting the conformal 
distance d{x\L) — L\ sin(7rx/L)|/7r, the leading term of the 
quartet correlations is of the bosonic form 



Q{x) =/9o 




cos(7rna;/2 + S) 
[d{2{m-x)\2L)]K/^ 

y/d{2x\2L)d{2m\2L) 
d{x + m\2L)d{x ~ m\2L) 



2/K 



(8) 



If one writes the Qi operator in a density-phase representation 
'^'), the first term is ^qiqi+x where qi — (QlQi) de- 



Vqie 

notes the local density of quartets (bosons). Because of OBC, 
Friedel oscillations appear close to the edge, leading to a typ- 
ical 2kp cosine term that decays algebraically /rom the edge 
(see a discussion in Sec. 13.6b . In terms of bosons, this decay 
of the density fluctuations 1 55 1 is controlled by Kb which gives 



o5 



< 
0.1 

0.01 


C: 


- ^=2.09 

- K=23l 

- K=236 
o DMRG 





8 
X 



32 



64 



Fig. 1. Typical fit of the quartet correlations in the MS phase on the 
SU(4) line. Data are obtained by DMRG with L = 128 and [//f = -4 
at filling n = 0.75. See text for the three fitting functions. Using 
Eq. ^ gives the Luttinger parameter K = 2.36. 
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Till 



Fig. 2. Density structure factor A*'(fc) obtained from QMC SSE sim- 
ulations at low temperature (T — 0.02t), fixed chemical potential 
H = — 2t and for various sizes L. Since finite-size effects are rather 
small, an accurate estimate of K can already be obtained on small 
systems from the small-fc linear behavior Here, for U/t = — 1 and a 
density close from 1.2, one gets K ~ 1.3 which is compatible with 
the DMRG estimate. 



K/4: for quartets. Note that the scaling dimension of the density 
operator is also K/A. Thus, our fitting procedure stems more 
from the phenomenology of hard-core bosons than from an 
exact result. Note that oscillations cannot be explained by the 
harmonics of the quartering correlations as derived in Eq. ^ 
because the exponents of the sub-leading terms are of order 2 
and 5, which are far too large to explain the oscillations (which 
actually increase with x). The amplitude co and phase-shift S 
are unknown parameters. The second term is the leading al- 
gebraic decaying term x^^^^ modified by finite size effects. 
The function replacing 1/x in between the brackets accounts 
for the vanishing of the wave-function at the edge of the box 
which induces a drop in the correlations. Typical data in the 
MS phase are given in Fig. [T] for a chain with size L — 128. 
Rather strong Friedel oscillations are observed in the signal (in 
the BCS phase, these oscillations are much smaller). Three fits 
are used to extract K. Firstly, a simple algebraic fit (which 
corresponds to taking cq = and L = oo in Eq. dHJ) yields 
K = 2.09. Secondly, a fit without the 2kp oscillations (cq = 0) 
gives K = 2.31. Thirdly, a fit using Eq. (|8]l with po, co, S and 
K as free parameters gives K = 2.36 and an excellent agree- 
ment with the data. It is thus important to take into account the 
finite size effects to have a reliable evaluation of K. In a pre- 
vious work 1 43 1, we have used an averaging of the correlators 
over i; this suppresses the oscillations but gives a less accurate 
estimate for K. A similar fit function as in Eq. (O but with a 
phenomenologically introduced cosine oscillations was used in 
Ref. 1 26 1 and leads to results very close to the ones obtained 
from Eq. (01. 

Another systematic way of calculating the Luttinger expo- 
nent is to use the density correlations N{x) and the associated 
structure factor N{k), where k is the wave-vector. This method 
is particularly suited for QMC as the density operator can be 
more easily sampled than the quartet operator. The value of K 
is extracted from the small wave-vector behavior of N{k): 



D ADW O BCS 
■ MS O MOW 




I O O ®-4 O O ■ O • 

Fig. 3. Phase diagram of the spin-3/2 Hubbard model ^2} at the in- 
commensurate filling n — 0.75, and for attractive interactions (U < 
or y < 0) from DMRG calculations (see text for definitions of the 
phases). The U — V line is the perturbative estimate for the transition 
between the two regions BCS/MDW and the ADW/MS regime. The 
dashed lines are the perturbative estimates for the crossovers between 
respectively BCS/MDW and ADW/MS. 



with a factor 4 in the denominator corresponding to the num- 
ber of fermionic flavors. For instance, this procedure has been 
shown to be very accurate for the spin- 1/2 Hubbard model [38l. 
An example of a typical fit is given on Fig.|2]at fixed chemical 
potential /i. Note that, since the QMC SSE algorithm is grand 
canonical, the density will slightly vary when the parameters 
(temperature or size) are changecfl. This effect can be seen from 
the position of the 2kp = TTn/2 peak in Fig.|2] One advantage 
is that a linear fit is simple to perform. However, in the limit 
of small densities, the 2kp peak approaches which makes it 
difficult to find the linear small-fc regime on finite size systems. 



3.2 Phase diagram at the generic density n = 0.75 

Fig. [3] displays the phase diagram of the spin-3/2 Hubbard 
model for attractive interactions (U < or F < 0) and a 
density n ~ 0.75. The density is chosen in such a way that 
the MS exists (from Ref. |26| we know that this is the case on 
the SU(4) line), and that there are no commensurate phased (in 
contrast to n = 1 1 15 , 43] or n = 2). Note that the MS phase is 
not accessible at filling n — 1 while the perturbative estimate 
of Eq. (|3]l predicts its existence |i43|. From its wide extension in 
Fig.|3] we observe that the MS phase, is very robust under the 
symmetry breaking term V. Thus, the quartet molecular phase 
is not an artifact of the SU(4) symmetry. This is an important 



4 fc^o fc ' 



(9) 



' Note that canonical algorithms are also available for such models. 

" Actually, one could argue that n = 3/4 is a simple fraction and 
that commensurate phases can occur. However, in terms of bosons, 
this would correspond to a density 3/16 for which a Luttinger ex- 
ponent Kb — 2/16^ is required to drive the transition |52|, giving 
K — 1/32 which is very small, but could, in principle, still appear at 
very large interactions. 
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result since in most of the realistic situations, the actual sym- 
metry is expected to be smaller than SU(4). Part of the answer 
is given in ID systems by the accepted view that, at sufficiently 
low energies and for generic interactions, the dynamical sym- 
metry is most likely to be enlarged [59 1: though the SU(4) sym- 
metry is not an exact symmetry, it is physically meaningful as 
an effective low-energy theory. As a consequence, the SU(4) 
model studied in Ref. [26J is a very good starting point to ex- 
plore the main features of the quartet phase. 

As a remark, we argue that the quartet phase also emerges 
in a problem with no extended SO(5) symmetry. For instance, 
we can consider the two-leg ladder model described in Sec. 12.21 
on the SU(4) line. The model simplifies to a model of two spin- 
1/2 Hubbard chains coupled only by the inter-chain density- 
density interaction V^. If we relax the constraint V\^ = U 
and let V± vary, we have the following picture. For Vj^ = 0, 
the SU(2) Hubbard chains are exactly solvable by the Bethe- 
ansatz technique and the corresponding Luttinger parameters 
"su(2)i ^su(2) can be computed exactly |60|. Bosonizing the 
V±_ coupling between the two chains gives, for the symmetric 
combination of the modes, the Luttinger parameter 



K+ = K 



SU(2) 



// 



1 + V^K, 



SU(2 



/7rwsu(2 



(10) 



K+ can be identified with K when Vj^ — U and governs the 
quartetting correlations as in Eq. (|5]l. We see that a negative 
Vj_ gives K+ > A'su(2)- Yet, we know from Ref. f60| that 
^su(2) ^ 2 in the limit of low densities and negative L/su(2) • 
A finite negative V± should thus easily stabilize a superfluid 
quartet phase with K^ > 2. Furthermore, we must note that 
Eq. ( [Tot is perturbative in V± but otherwise valid for arbi- 
trary values of usu(2) ^nd Ksij(^2)^ which are known even in 
the strong coupling regime |C/su(2)l/^ ^ 1- In particular, the 
fact that K saturates along the U/t line of the SU(4) Hub- 
bard model |26| might be similar to the saturation of Ksu(^2) 
at large |C/su(2)|- For both models, the saturation is associated 
with the onset of a hard-core boson regime of pairs for SU(2), 
and quartets for SU(4). A similar superfluid phase where pair- 
ing of bosons along the rungs occurs in the bosonic ladder 
model f6r|, the phenomenology and perturbative argument be- 
ing essentially the same. 



3.3 Excitations gaps and spectra 

We now turn to the excitation gaps in the two regions of the 
phase diagram. In the context of cold atoms experiments, the 
quartet phase can be probed by radio-frequency spectroscopy 
to measure the excitation gaps of the successive quartet disso- 
ciation process. In particular, the existence of molecules can be 
characterized by finite one and two particle gaps while four 
particle excitations remain gapless. This feature alone does 
not distinguish between various phases (dominant superfluid 
or density correlations) but it allows to check for the formation 
of bound-states. 

The energy gap to fill when adding p particles in the system 
is defined by 



O f//(=-l 
□ f//(=-2 
A U/l=-4 



P^ 



4p - 
0.1 



ff 







z\. 



EoiNf +p) + EoiNf+p) - 2EoiNf) 



(11) 



0.01 0.02 0.03 

1/L 

Fig. 4. One, two and four particle gaps along the SU(4) line for the 
density n = 1. Left: scaling of the four particles gap for different U. 
Middle: Aip and A2p as a function of U/t on a finite system with 
L = 12, and extrapolated Z\4p. Right: We also show the inverse of 
the correlation lengths 5ip,2p obtained from the Green's function and 
the pairing correlations. Inset: comparison of the numerically obtained 
ratios and the y/2 prediction. 



with EQ{Nf) the energy of the ground-state with Nf particles. 
We choose Nf = 4(2to +1), with m an integer, so that we 
would have closed shells in the case of periodic boundary con- 
ditions. Fig. |4] provides the results on the one, two and four 
particles gaps on the SU(4) line with a density n — 1 and a 
small system L = 12. For large sizes and interaction \U\, quar- 
tets are strongly bound and convergence of systems with Nf 
not a multiple of four can fail. Indeed, even after many sweeps, 
the density distributions are not symmetrical with respect to 
the center of the chain. Results given in Fig. |4] are those with 
symmetrical ground-states and well-converged energies. The 
small size of the chain may cause finite size effects but, still, 
a clear opening of the one and two particle gaps is found. For 
A4p, no convergence issues are found and scaling can be per- 
formed. Fig. m shows that A4J, = in the thermodynamical 
limit, in agreement with the algebraic decay of the correlations. 
To further check the consistency of the numerics and the low- 
energy theory, we use the fact that the ratio Z\2p/Z\ip is known 
to be exactly \/2 from the integrability of the SU(4) Thirring 
model 1 62 1 describing the spin part of the Hamiltonian in the 
low-energy approach. This ratio is plotted for the L — 12 gaps 
in the inset of Fig. |4] and agrees reasonably well with the pre- 
diction. 

Another way of probing the presence of finite gaps is to 
look at the associated correlation functions. The one and two 
particles gaps, when finite, are associated with the exponen- 
tial decay of the Green's and pairing correlation functions, as 
observed in Ref. [43 1. In this case, the two coiTelation lengths 
behave as ^ip,2p ~ M/^ip.2p, where u is the sound velocity of 
the bosonic mode. The inverse correlation lengths are given as 
a function of interaction in Fig. H) The same universal ratio is 
expected for the correlation length and well observed numer- 
ically on Fig. m These results support the correctness of the 
low-energy approach, even for strong couplings, and in partic- 
ular the validity of the separation of the spin-charge sectors. 

Similarly, gaps and the inverse correlation length can be 
computed by DMRG on the BCS Une. Results are given in 
Fig. |5] which shows a smooth opening as soon as interactions 



G. Roux et al. : Spin 3/2 fermions with attractive interactions in one dimension 




Fig. 5. One and two particle gaps along the BCS line for the density 
n = 1. Left: Scaling of A2p. Middle: Aip on a finite system with 
L = 12, and extrapolated A2p. Right: the inverse correlation length 
of the Green's function ^ip is also shown. 




Fig. 7. Fraction of molecules in the system as a function of U for the 
SU(2) and SU(4) models as defined by Eq. l ll2b from a system with 
L = 64 and density n = 1. 




Fig. 6. Lowest energy E{k) as a function of the momentum k for a 
ring of length L = 8 with 8 fermions for various interactions (a) along 
the SU(4) fine (V" = 0); (b) along the BCS line ((/ = 0). Antiperiodic 
boundary conditions are chosen to have closed shells. Energies are 
measured relative to the ground-state energy. 



are turned on. Two and four particle excitations are gapless 
while a one-particle gap opens. We only show Z\2p because 
Z\4p also scales to zero when Z\2p = in the thermodynamical 
limit. These behaviors follow the results obtained in Ref. [43! 
for the pairing and correlation functions. To give additional in- 
sights on the opening of the one-particle gap, we provide the 
evolution of the inverse correlation length of the Green's func- 
tion ^ip. Experimentally, the one-particle correlation length ^ip 
will appear in the momentum distribution of the condensate. 

Exact diagonalization on small systems on a ring allows 
for the computation of the excitation energy spectrum E{k) 
vs. momentum k. Even if finite size effects can be important, 
some qualitative information can be extracted. Figure |6] dis- 
plays the spectra along the SU(4) and BCS lines for the density 
n — 1. Note that anti -periodic boundary conditions are used to 
have closed shells when U = V — for a chain with 4(2m) 
fermions with m an integer On the SU(4) line, we observe that 
the 2kp excitation has a lower energy than the Akp, which is 
associated with the dominant density fluctuations at 2/ci? in this 
region of parameters. As |t/| increases, all energies go down, 
so the sound velocity of the charge mode u also decreases, in 
agreement with the perturbative estimate of Eq. (|4|i. The spec- 



trum evolves continuously towards the strong-coupling limit. 
On the contrary, along the BCS line, a crossover is found be- 
tween a regime, at low \V\, in which the minimum is at 2fci?, 
and the strong coupling regime for which the minimum is at 
4kp. We will see hereafter that a similar crossover is found 
in the density fluctuations and Friedel oscillations. Lastly, one 
can note that the sound velocity u slowly decreases through 
this crossover line (and slower than on the SU(4) line), again in 
agreement with the perturbative prediction. 



3.4 Quartets formation on the SU(4) line 

In this section, we discuss the crossover from the weak- 
coupling regime to the strong-coupling regime on the SU(4) 
line as the attractive interaction is increased. When |f7| is large, 
the physics is essentially the one of hard-core bosons with re- 
pulsive interactions, as it has been discussed in Sec. 12. II within 
the low-energy approach. To investigate how the quartets form, 
we can compute the local density of these "molecules". From 
a more general point of view, and to compare with the SU(2) 
case, we consider iV-particle bound states in the context of the 
SU(-/V) Hubbard model ll26l . The local density of molecules 
is 'm{x) = {nx,i ■ ■ -rix^N) (which we denote by q{x) for 
quartets). For free fermions, this operator has a finite expec- 
tation value that we subtract to keep only the connected part 
m{x) = {ux^i ■ ■ ■ tIx.n) — [n/N)^ . If molecules are tightly 
bound on-site, we expect {rtx^i ■ ■ ■ n^.N) to be close to n/N 
though slightly lower. Therefore, we can define a molecule 
fraction (number between zero and one) as 



N 



%Molecules = 



m{x) - {n/N) 
n/N - {n/N)N 



(12) 



where the bar means averaging over all sites. The evolution of 
this quantity along the SU(2) and SU(4) lines are compared in 
Fig.|2]as a function of NU (and not U) because the interaction 
term scales like N"^ while the kinetic term only scales as N . 
At large \U\, molecules are tightly bound, but the SU(4) model 
is closer to a hard-core boson model than the SU(2) one. An- 
other difference is the low-C/ increase which is linear for SU(2) 
and power-law for SU(4) with an exponent larger than 2. Note 
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Fig. 8. (Color online) Luttinger parameter K vs filling n for various 
U obtained by fitting quartet correlations as in Fig. [T] for the SU(4) 
model. K ^> A when n -^ and K ^r 1 when n -^ 2 are the 
asymptotic behavior of the strong-coupling limit. 



that the behavior should also depend on density, particularly at 
small U . 

Consequently, we expect the large negative U physics to be 
essentially the one of hard-core boson. Still, we emphasize a 
major difference between SU(2) and SU(4): from perturbation 
theory, the SU(2) case leads to hard-core bosons with equal 
effective hopping and nearest-neighbor repulsion (also equiva- 
lent to an effective spin- 1/2 XXZ chain ll63l ): on the contrary, 
in the SU(iV) case (with N > 2), the effective hopping at iV* 
order in perturbation theory behaves as t'^ /\U\^~^, so it is 
negligible compared to nearest-neighbor repulsion, which is of 
order iV|C/|. 



3.5 Evolution of the Luttinger parameter and the 
commensurate phase ADW for n = 2 

The previous considerations allow for a simple interpretation 
of the behavior of the Luttinger parameter X as a function of 
the density n for large negative U . DMRG results are given 
in Fig. [8] As expected from the strong coupling argument, K 
decreases from K ^ N io K ~ N/A (from 4 to 1) as the 
density n varies from to half-filling (n = 2). Particle-hole 
symmetry would give the behavior for 2 < n < 4. Again, 
we recall that molecular superfluidity is the dominant instabil- 
ity when K > 2, which is generically the case at low enough 
density. When n = 2, a fully gapped phase is obtained with 
short-range quartet correlations. The phase is two-fold degen- 
erate with a IT ordering of the local density (one quartet every 
two sites). Hence, we call this phase ADW^. In terms of an ef- 
fective bosonic model discussed in the previous paragraph, this 
corresponds to a "charge" density wave phase of the equivalent 
bosonic model at half-filling |64,65|. The density of bosons 
being 1/2, the corresponding critical value for their Luttinger 
parameter fS^l is Ki, = 1/2 (at fixed density, changing interac- 
tions), which gives K = 2 for our model. As K = 1 < 2, the 
ADW^ phase emerges as soon as the interaction U is turned on. 
Working at fixed interactions and varying the density, the criti- 
cal value is now Ki, = 1/4, which gives the observed limiting 
value K ~ 1 as one approaches n — 2 (see also Ref. [j26|). 



3.6 Density fluctuations and Friedel oscillations 

The density fluctuations can be analyzed from the correlations 
structure factor N{k) with QMC, and from the behavior of the 
Friedel oscillations of the local density in open chains, as usu- 
ally done in DMRG. The Friedel oscillations are the response 
of the fermionic density to the open end of the chain, which 
acts as an impurity. These modulations can give access to Lut- 
tinger parameters |66|. Data for the SU(4) line (not shown), 
and more generally in the ADW/MS region of the phase dia- 
gram, are all consistent with the low-energy predictions 1141 
[T5ll43l N{x) ~ cos{2kFx)x~^^^ for the density correlations 
and n{x) ~ cos{2kFx)x^^''^ for the Friedel oscillations. Note 
that the N{2kp) peak diverges with the system size L provided 
K < 2, signaling the quasi-ordering of the density fluctuations 
oftheADWphase|43|. 

Friedel oscillations in the BCS phase have a different be- 
havior. As shown in Fig. |9] for the generic density n = 0.75, 
there is a qualitative change in the wave-vector of the oscil- 
lations from 2k F at low |1^| to Akp at large \V\. The predic- 
tions are that the 2k p term should be short-range but a Akp 
term can develop with coiTelations N{x) ^ cos{4:kpx)x~'^^ . 
For the Friedel oscillations, we thus expect a leading contri- 
bution behaving as n{x) ~ cos{4:kpx)x^'^ , similar to what 
was found in two-leg ladders |67|. To explain the behavior ob- 
served in Fig.|9] we argue that at low \V\, the amplitude of the 
2kp term remains significant (it is finite for a free system at 
V ^ 0) and with a coiTelation length which is still large (see 
Fig. |5]l. When \V\ increases, the Akp term emerges with an 
increasing amplitude. Fits have been carried out in Fig. |9] us- 
ing n{x) = no + ni cos{2kpx + 6)e^^/^ for V/t — —0.5 
and n{x) ~ uq + niCos{4:kpx + S)/[d{x\L)]^ for larger 
\V\. The Luttinger exponents obtained from the fits are close 
to the ones obtained from the pairing coiTelations. In addition 
to DMRG calculations, DQMC data (see for instance Fig. 4 of 
Ref. P3l and other results not shown) support a similar quaU- 
tative change in the wave-vector and no divergence of the Akp 
amplitude with the system size. Indeed, this divergence only 
occurs for K < 1/2, i.e. in the MDW phase. Lastly, the same 
crossover around V/t = — 1 is found in Fig.|6jb). 

The crossover to the large \V\ physics can be qualita- 
tively understood within the following picture: when V is large. 
Cooper pairs have a tendency to form on-site, and certainly 
repel each other to gain local kinetic energy. This gives a 
typical Akp fluctuation of the local density and kinetic en- 
ergy as found in Fig. |9] The local kinetic energy term is 

ii^) = (Z]cr4+i<TCxCT>. In the ADW/MS region, it follows 
the Friedel oscillations of the density, so we have t{x) ^ 
co8{2kpx)x-^/'^. In the BCS/MDW region, the 2kp compo- 
nent is short-range so the leading term is the Akp one, t{x) ^ 
cos{4:kpx)x^^ , as for n{x). The enhancement of these fluc- 
tuations (the total kinetic energy rather decreases with interac- 
tions as seen in Fig.|9]l as \V\ is increased is reflected through 
the decrease of the Luttinger exponent K, as it was found for 
n = 1 in Fig. 5 of Ref. |43|. A similar slow decrease with val- 
ues of ii' lower than one at large \V\ is found for n = 0.75. This 
decrease is not predicted in the perturbative estimate of Eq. ^ 
and is therefore a typical strong-coupling behavior. If one adds 
the repulsive interaction U on-site, strictly on-site pairs are no 
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Fig. 9. Local density n{x) and kinetic energy t{x) along the BCS line 
at filling n = 0.75 on a L — 128 chain for increasing \V\. Friedel 
oscillations of the density have been shifted vertically for clarity. Thin 
lines are fits to results from DMRG calculations (see text). 



more favored and the pairs lower their energy by delocalizing 
themselves on a bond. If U is large enough and the density 
commensurate at n ~ 1, this qualitative picture leads to the 
bond-order wave phase observed in Fig. 7 of Ref. i43l which 
breaks translational symmetry and is two-fold degenerate. 



3.7 The transition between BCS/IUIDW and ADW/IUIS 

This section gives some results on the transition line between 
the two regions BCS/MDW and ADW/MS. It was already 
shown that it belongs to the Ising universality class and that 
the ratio between the pairing and quartet correlations R{x) = 
P'^{x)/Q{x) has the universal behavior 1/x at the critical 
point in agreement with conformal field theory (CFT) predic- 
tions lfT5]lfT]l43ll. 



3.7.1 Effect of the density on the transition line 

We first investigate the question of the dependence of the tran- 
sition line with respect to the density n. This is an important 
issue for inhomogeneous systems, such as trapped cold atoms, 
since the local density in the cloud evolves continuously from 
zero to a finite value in the bulk. We know that the perturba- 
tive result U = V for this line does not depend on n, while the 
crossover lines in Fig.[3]noticeably depend on n from Eq. ^. 
In the strong coupling regime, we study numerically the transi- 
tion line for several densities. To that purpose, we first use the 
behavior of R{x) averaged over distances ranging from x = 20 
to 40 (from correlations in a system with L = 128), which we 
denote by R. In Ref. |43|, the same ratio was used for a fixed 
distance x — A5 for n — 1; the two procedures lead to the 
same results but the first one is more suited to systems with 
low densities as oscillations have a longer wave-length. From 
Fig. [TOl' a). we see that R vanishes linearly around the critical 
point which is due to the Ising nature of the transition [43 1. 
However, the higher the density, the wider the range of U over 
which the linear behavior is observed. At density n = 0.25, 
the linear behavior is not recovered, certainly because of our 




-1.2 -1 

U/t 

Fig. 10. Effect of the density on the transition line with fixed V/t = 
—2 (a) Averaged ratio R signaling the BCS/MS transition for different 
densities. The critical point U/t ~ —1.2 of the BCS/MS transition 
hardly depends on density, even in the strong coupling limit, (b) The 
fraction of molecules (quartets) does not give a direct estimate of the 
transition but its derivative (inset) does. Open symbols are for a system 
with L — 128 while symbols filled with grey correspond to i = 32 
(finite size effects are stronger at large density but remain small). 



mesh points. The main result is that the position of the criti- 
cal point U/t ~ —1.2 hardly depends on the density, though 
corresponding to the strong coupling regime. 

Experimentally, it would be impossible to measure R so it 
is interesting to look at the behavior of the molecule fraction 
through the transition. This is given in Fig. fTOl b) in which one 
finds that the molecule fraction increases from the BCS/MDW 
to the ADW/Ms region, but rather smoothly at large densities. 
However, taking the derivative of the curve with respect to U 
gives a clear peak pointing at the critical point. Consequently, 
this would be a rather simple way to locate the transition line 
experimentally as it could be measured and works even for an 
inhomogeneous cloud. The dependence of the peak with re- 
spect to the density is found to be small. 



3.7.2 Entanglement entropy and central charge 

Another prediction from CFT concerns the behavior of the cen- 
tral charge c of the model at the transition line. The central 
charge somehow measures the effective Ising degrees of free- 
dom in the low-energy physics. It is expected to be one in both 
regions (one gapless bosonic modes) around the transition but 
exactly equal to 3/2 at the transition due to the emergence of 
the Ising criticality with central charge c = 1/2. A simple way 
to extract the central charge with DMRG is to use the behavior 
of the von Neumann entanglement entropy Syt.i{x) of a block 
of size x < L.lt is defined as 



Syti{x) ^ -TT[p{x)lnp{x)] 



(13) 



where p{x) is the reduced density matrix of the block. As 
has been emphasized recently in several studies, the use of 
the entanglement entropy can provide crucial information for 
condensed matter studies since it allows to detect quantum 
phase transitions without any knowledge on the order param- 
eters [68]. It is straightforwardly computed with the DMRG 
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Fig. 11. Von Neumann block entropy Syti{x) for a block of size x 
and local kinetic energy t{x) around the critical point U /t — —1.2 
for fixed V/t = -2 and n = 0.75 (see Fig. QHll with L = 128. 
Fits (filled circles) using Eq. ( I14l l are quite accurate, allowing for the 
determination of the central charge c. 
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Fig. 12. SU(4) model: N{k) obtained from QMC SSE simulations 
at various temperatures (in units of t) for U /t — —2, L — 32 and 
jj, = — 3.4f . At low temperature, the small-fc linear behavior allows to 
extract K ~ 1.6 which is compatible with DMRG estimate of 1.5 for 
same parameters. Moreover, the temperature effect allows to estimate 
at which energy scale the 2kF peak appears, kp shown on the plot 
corresponds to the mean density (n) = 1.4. 



algorithm from the eigenvalues of the reduced density matrix 
which is obtained at each iteration. Following the ideas devel- 
oped in Refs. II691I70L a subleading oscillating term emerges 
due to open boundary conditions. Similarly to what happens 
in the XXZ model, we expect the oscillations to be related to 
the local kinetic energy t{x) (similar to the dimerization term). 
Indeed, we can carry out fits in both regions by using the fol- 
lowing ansatz 



S,j^{x) = - In d{x\L) + A + B{t{x) - t) 



(14) 



where t is the mean value of t{x) in the bulk, and A, B are two 
constants. The behavior of the entanglement entropy thus gives 
access to the central charge c. On Fig.[TT] a clear jump of c is 
observed at the transition. The expected value c = 1 is well 
reproduced in the two regions away from the critical point. At 
the transition, the fit is not as accurate but restricting it to the 
bulk region yields c == 1.51, at the price of describing less 
well the strong oscillations close to the edges (a fit including 
all data yields c = 1.71 but overestimates the behavior in the 
bulk). Note that the singular behavior of c differs from the con- 
tinuous behavior of K and the gaps at the transition. Note also 
that the ikp oscillations in the ADW/MS region are reminis- 
cent of similar features seen in local density and kinetic energy, 
and compatible with a 2/cf soft mode (see Fig.|6h) |71 1. A last 
remark is that the local kinetic energy oscillations on the tran- 
sition line should follow t{x) ^ coa{2kFx)x^''-^^^'>^'^, i.e. an 
exponent between those in the neighboring BCS/MDW region 
(provided K > 1/3) and ADW/MS region. 



4 Experimental signatures of MS phase 

This section is devoted to the study of effects particularly rel- 
evant for experimental set-ups. In addition to these results, we 
have already discussed in Sec. l3.3l the excitation gaps, relevant 



for radio frequency spectroscopy, and the molecule fraction, 
which could be measured experimentally from pictures resolv- 
ing the hyperfine states. 



4.1 Effect of temperature on the 2kF peak 

Using QMC, it is possible to investigate the energy scale 
at which the zero-temperature features become relevant. In 
Fig-El we study the effect of temperature on the density struc- 
ture factor N{k) calculated on the SU(4) line for different tem- 
peratures. Let us remind that we work in a grand-canonical en- 
semble so that the density (n) varies with temperature (typi- 
cally between 1.2 and 1.6 for this plot). In particular, we have 
shown the 2k p location corresponding to the low- temperature 
density ((n) = 1.4). The two properties of interest are the 
emergence of the 2kp peak, and the linear behavior at small- 
k. We observe that, below a typical temperature of order O.lt, 
these two features qualitatively approach their T = behavior, 
while a quantitative estimate requires a much lower tempera- 
ture of order 0.02t. If these features could be measured exper- 
imentally, one could identify the two main regions from the 
strength of the 2kp and Akp peaks. 



4.2 Effect of the trap 

The effect of the trap confinement on two-component 
fermionic gases loaded in optical lattices has been studied 
for repulsive [72] and attractive II63II73I I interactions. Trapped 
fermions with attractive interactions were also studied for im- 
balanced populations ll74l . i.e. with no SU(2) symmetry. The 
Hamiltonian term corresponding to the harmonic confinement 



^^(*-(i + l)/2)V, 



(15) 
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Fig. 13. Typical density profiles of the trapped condensate for fixed 
"density" Nfui = 0.2 and U/t = —4, with respectively Nf — 
36, 68, 100, 140, 180. The inset shows the scaling of the Wlo in the 
thermodynamical limit, which well agrees with the TG behavior. 



with uj the trap frequency and {L+l)/2 the middle of the chain. 
We thus take the box width L to be larger than the cloud's width 
not to induce boundary effect from the edges of the box. With- 
out a lattice, the chemical potential reads /i ~ ujNf/N for free 
fermions. The thermodynamical limit is understood as taking 
the w ^ limit while keeping Nfuj constant so that the den- 
sity at the center of the trap remains constant. Consequently, 
NfU! is similar to an effective density of the system and, de- 
pending on it, several regimes are identified. 

The density profile of the condensate n{x) (or rii — n{xi) 
in case of a lattice) is directly accessible from experimental pic- 
tures. For free fermions without an optical lattice in the Tonks- 
Girardeau (TG) regime with iV-color I175II76II771I . one has 



njcix) == rioA/l - a;2/i?|, 



(16) 



with Rjc ^ y/Nf/Nu! and uq ^ y/NfUi/N, on which 
we observe that uq is kept constant in the thermodynami- 
cal limit. These results are valid if the trap evolves smoothly 
enough such that the local density approximation (LDA) is ex- 
pected to be a reasonable assumption. In presence of an op- 
tical lattice, the energy per particle and the dispersion rela- 
tions are changed. In this situation, LDA gives a density pro- 
file of the type n{x) = no arccos(x^/i?^ — b'^)/ arccos(— fe^) 
for 



I a; I < i?Vl + 62 fJW\. The typical width of the den 
sity distribution, which can be measured, is defined as M^ = 

2 



Nf 



^j(i — io)2rii. In the Tonks-Girardeau regime, we 

have W ^ Rjg "^ ^/NfJNuj, so that W ^ uj^^ in the thermo- 
dynamical limit. We found a similar behavior for the spin-3/2 
fermion model under study with a scaling which agrees well 
with the TG one, as one can infer from the results of Fig. [13] 



4.2.1 Atomic density waves 

For sufficiently smooth traps, the bulk of the condensate fea- 
tures the typical 2kp oscillations reminiscent of the ADW 



T 



n ' r 



n ' r 




Fig. 14. Typical density profile of the trapped condensate on the SU(4) 
line (U/t = -4) with a smooth trap (Nfu = 0.2 and Nf = 140). 
One observes that 2kF — 7r/2 oscillations develop at the center of 
the cloud. Close to the edges, the wave-length of the oscillations in- 
creases, qualitatively following the local decrease of the local aver- 
aged density. A reasonable fit for the bulk physics is obtained using 
the Tonks-Girardeau result plus an oscillating term. 



phase encountered with open-boundary conditions. For in- 
stance, Fig.[T4lshows a typical density profile on the SU(4) line. 
The density profile can be reasonably fitted up to the edges by 
a Tonks-Girardeau profile (Eq. ( fTSI l) plus an oscillating term 

n{x) — njcix) + Sn cos{2kp {x)x) , (17) 

with the effective Fermi wave-vector 



kpix) = -jnoyl - x^/{Rtg) 



(18) 



In Fig. [141 we first fit the TG profile and subtract it from the 
data to only keep the oscillating term that we fit using the same 
value of riQ. The slight dependence of kp on x accounts for the 
increase of the wave-length as the density decreases towards 
the edges of the condensate. However, if Rjg is a free param- 
eter of the fit, we find that Rjc — 2i?TG gives a better fit of 
the oscillations. This discrepancy could be due to finite size 
effects, as the TG profile should be valid for large enough sys- 
tems and far enough from the edges of the condensate. The 
condensate has sharp edges (LDA usually fails to explain the 
behavior close to the edges) and in the following, we call a the 
radius at which the density vanishes (a < i?TG)- 

The main question is now to discuss the thermodynami- 
cal limit of the oscillations amplitude in the bulk 6n. For non- 
trapped gases in a box, we expect the oscillations to be zero in 
the middle of the system (except for a translationally breaking 
phase) as the Friedel oscillations decay away from the bound- 
aries. In the SU(2) case, a finite Sn has been found around the 
commensurate density uq — 1 063II73I . These atomic-density 
waves will have clear signatures in the density structure fac- 
tor that can be measured with light-scattering diffraction. The 
latter is defined by 



S{k) 






(19) 



12 



G. Roux et al: Spin 3/2 fermions with attractive interactions in one dimension 




Fig. 15. The density structure factor S{k) corresponding to Fig. [T3l 
(same color code), and the scaling of the 2kF peak position and am- 
plitude in the thermodynamical limit. 
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Fig. 16. Correlations in a trap system with u) = 0.002, f/ = — 4, V = 
and Nf — 36. (a) Quartet correlations Q{x) showing dominant 
MS correlations (from comparison with 1/x) (b) Pairing and Green's 
function are short-range. The full (open) circles account for the +(-) 
sign coming from the sin(fcFa;) term. 



and an example is given in Fig. [15] The main features of the 
structure factor is the peak at 2kp ~ 7rno/2 which signals the 
ADW oscillations of the trapped systems. Oscillations in S{k) 
are due to the finite size 2a of the condensate and may van- 
ish for large systems. We expect the height of the peak to be 
proportional to ((5n)^ for large enough systems. The insets of 
Fig. [15] show the evolution of ^S'(2fci?) "^ Sn and 2kp (de- 
fined as the k at which the peak has its maximum) as a function 
of Nf . Fits are obtained for both quantities with an exponential 
law So + sie'^f^^. We infer from these results that 6n is fi- 
nite in the thermodynamical limit for our choice of parameters, 
similarly to what was found in the SU(2) case 



4.2.2 Correlations and the BCS/ADW transition 

We now turn to the behavior of the correlations in trapped 
gases. They are computed by fixing one point of the correla- 
tor at the center of the cloud and letting the position x of the 
other ranging from the middle to the edge of the condensate. 
For a trapped condensate, a gap computed from energy differ- 
ences can be spoiled by edges effect. To check the gapped na- 
ture of the one- and two-particle excitations along the SU(4) 
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Fig. 17. (a) Correlations in a trapped system with U = —4 and V — 
(a) for fixed lu = 0.002 and various Nf. The distance x is rescaled by 
the trap radius a. (b) The same for fixed Nf = 36 and varying id. 



line, we compute the Green's function and the pairing corre- 
lation. From the behavior of the quartet correlations, we can 
deduce the leading fluctuations in the bulk of the cloud. These 
results are given in Fig.[T6]in which we see that for this small 
value of the effective density (N fuj — 0.07), the physics in 
the bulk is essentially the same as for the non-trapped sys- 
tem. Quartet correlations slightly increase as one approaches 
the edge of the condensate. Assuming the LDA approxima- 
tion to work in such small systems, we expect the local Lut- 
tinger parameter M79II80I controlling the power-law decrease of 
the correlations K{x) = K{n{x)) to increase with decreasing 
densities, following the homogeneous system result of Fig. [8] 
As the quartet correlations behave as 2/ K and because K de- 
creases with density, they are actually reinforced by the har- 
monic confinement. Trapping thus favors the observation of the 
MS phase. By lowering the effective density NfU, we expect 
a crossover from sub-leading to leading quartet correlations. In 
Fig. [17] we show both the effect of changing the number of 
particles Nf while keeping lo constant and the opposite situ- 
ation where lo varies. We find that quartet correlations always 
decreases slower than 1/x (corresponding to the critical value 
Kc = 2) below approximately NfUj ~ 0.08. 



4.2.3 Effect of varying interactions and deep trap physics 

In this section, we address the question of the effect of varying 
interactions on the density profile of the condensate. We expect 
the width of the condensate to strongly depend on interactions: 
repulsive interactions make the condensate inflate while attrac- 
tive interactions can strongly reduce W. In the large trap fre- 
quency limit, we furthermore have a minimal width of the con- 
densate due to Pauli's exclusion principle which depends on the 

number of colors: VT™" = J\{{Nf/NY- l) - Nf/y^N. 
Thus, starting from free electrons and keeping NfUj/N con- 
stant, the ratio W^'^/W^""" ^ y/Nfbj/N should be constant. 
We give a comparison with the evolution of the width in the 
SU(2) case as a function of NU: the collapse of the conden- 
sate is faster in the SU(4) case. In Fig. [18] we show, for a con- 
stant and rather large effective density Nfuo — 1.4, the evo- 
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Fig. 18. (a) Typical density profiles of the trapped condensate for var- 
ious interaction strength U on the SU(4) line. There are Nf — 28 
fermions and uj — 0.05. (b) Width of the condensate as a function of 
U/t. Its width is roughly divided by a factor 2.5. The SU(2) line has 
been computed with the same u but with 14 fermions. The width is 
much more sensitive to U in the SU(4) case, (c) Structure factors have 
been rescaled for clarity. 
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Fig. 19. Evolution of the profiles along the V — —2 line with decreas- 
ing U showing the crossing of the transition (located at U/t = —1.2). 



lution of the condensate when increasing (negative) U along 
the SU(4) line. Two effects are observed: the width of the con- 
densate strongly decreases which induces a sharp increase of 
the density at the center of the trap. Consequently, the effec- 
tive Fermi wave-vector increases, which can be observed in the 
density structure factor. Similar effects have been found for the 
SU(2) attractive model 1 8 1 , 82 1 . Because the effective density is 
large enough, the commensurate ADW phase develops when 
no > 2. A strong signal at fc = tt is then present in the density 
structure factor. 

Moving away from the SU(4) line, we can investigate the 
behavior of the density profile across the transition between the 
two regions of the phase diagram of Fig. |3] Indeed, we have 
seen in Sec. 13.7. fl that the transition line hardly depends on the 
density. Therefore, for inhomogeneous density profiles, we ex- 
pect the transition to hold for identical parameters, and that the 
whole cloud will be affected by the phase transition (meaning 
that no domains will form). Note that this is not true for the 
crossover lines inside each region of Fig. [3] and that is seen 
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Fig. 20. Starting with a system with strongly bound quartets (TV/ = 
28, U/t = —8), one can increase the frequency of the trap. Quartets 
have a tendency to repel each other but when the trap is too deep, they 
progressively melt at the center of the trap. 



through the variation of K{n{x)) in Fig.[T7l The typical evo- 
lution ofn{x) along the V/t = —2 line (varying U) is given in 
Figim Two main effects are visible: the width of the conden- 
sate is roughly divided by two, and the amplitude of the Friedel 
like oscillations are strongly enhanced in the bulk. Remarkably, 
because of the strong increase of the density at the center of the 
cloud in the ADW/MS region, the wave-vector of the oscilla- 
tions hardly changes (it is ikp in the BCS phase, but 2kp in the 
ADW/MS region). If one is able to measure the local density of 
molecules, a crossover from zero to a significant value will be 
also observed across the transition, similarly to what was found 
in Fig. [TOl The fact that the density of quartets is spread over 
the whole condensate through the transition line supports the 
absence of domain formation. 

If the trap is very deep, corresponding to large values of 
the effective density NfOj, all atoms will form an homogeneous 
condensate in the middle of the trap, the width of which is H/™'" 
as discussed above. As displayed in Fig.|20] this state emerges 
from the melting of the ADW" phase which can be qualita- 
tively understood as the mere competition between the effec- 
tive nearest neighbor repulsion of the quartets and the potential 
energy of the trap. This effect is very similar to the one found 
for the SU(2) model in Re f. 1821 . Forui — 0.05, the density pro- 
file is shifted from the center of the trap. Actually, the energy 
of such a shifted state is much smaller than other energy scales 
(equal to 0.00875t within a classical approximation), so that it 
is nearly degenerate to the ground-state and DMRG gets locked 
into it because the effective hopping term of the molecules is 
too small (for large \U\). 



4.3 Flux quantization 

In standard electronic systems, flux quantization experiments 
can directly measure the electric charge of the carriers by con- 
sidering a ring geometry threaded by a magnetic flux and look 
at the flux periodicity of the total energy. In the case of neu- 
tral cold atoms, such a flux analogy could be realized thanks 
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Fig. 21. Energy vs flux for a periodic chain of length L = 8 with 8 
particles for various interactions. Energies are measured relative to the 
absolute ground-state energy (here at <^ = tt). (a) SU(4) model show- 
ing the appearance of 4 minima; (b) BCS phase with 2 minima. Note 
that only flux between and tt are shown since data are symmetric 
E{—(j>) = E{4>) and 2it periodic. 



to the possibility of rotating the trap. Therefore, it could be 
possible to prove the existence of N particle bound states by 
checking if minima of the energy are degenerate |T5l. More- 
over, ring-shape geometries have been realized experimentally 
Il83ll84|[85l so that such experiments could be performed in the 
near future. 

We have performed exact diagonalization for the SU(4) 
case on small chains of length L = 8. Although these sizes 
are relatively small, we expect that bound-state formation can 
already be checked since it is a local process. Indeed, as the 
interaction strength increases, we see that the number of min- 
ima changes from one to four (see Fig. l2lT a)) in the ADW/MS 
phase, while it increases from one to two in the BCS/MDW 
phase (see Fig. lSTT b)). in full agreement with the predictions of 
the low-energy approach ifTSl . These observations are compati- 
ble with our predictions of four- and two-particle bound states, 
respectively. Note that the overall energy scales decrease since 
the band width is proportional to the effective molecular or pair 
hopping. Similar calculations suggesting pairing and superflu- 
idity in a multicomponent fermionic model have been proposed 
recently [86 1. As a final comment, let us remind the reader that 
the effect of such a flux corresponds to a twist in the boundary 
conditions, and therefore vanishes in the thermodynamic limit 
where E{(f)) becomes flat. Still, such an effect can be observed 
on finite lattices. 



5 Conclusion 

Motivated by fermionic cold atoms experiments where generi- 
cally many hyperfine states coexist, we have investigated spin 
3/2 fermions with contact interactions in an optical lattice. We 
focus on the attractive case for a generic density. By using 
large-scale numerical techniques, we describe the phase dia- 
gram and discuss all competing phases. In particular, at low 
density, we confirm the existence of a large molecular su- 
perfluid phase where dominant correlations are superfluid-like 
made of four-particle bound-states. This phase has a large ex- 
tension and is not restricted to the SU(4) model. In another re- 
gion of the phase diagram, two-particle pairs become gapless. 



leading to a BCS phase. We have shown that the phase transi- 
tion between MS and BCS phases can be located using entan- 
glement measurements, such as the von Neumann entropy or 
the molecule fraction. 

In order to make contact with possible experimental obser- 
vations of such phases, we have investigated the role of the 
trapping potential. In many respects, correlations inside the 
bulk of the condensate are similar to the homogeneous case 
if the effective density is low enough. ADW oscillations can be 
probed from the density structure factor. Furthermore, we give 
an estimate for the crossover effective density below which 
leading MS fluctuations are dominant. Moreover, playing with 
the trap can provide useful informations about the size of the 
condensate and the density correlations, which are accessible 
experimentally. For instance, deep in the ADW/MS phase, the 
physics can be understood from tightly bound "molecules" that 
act as hardcore bosons. These objects could be measured ei- 
ther by looking at the molecules fraction or by using rf spec- 
troscopy. Finally, we propose to distinguish between MS and 
BCS phase by using the molecule fraction or ring-shape ge- 
ometries. We hope that such experiments will be performed in 
the near future. 
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